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Abstract 

In this paper we re-examine the traditional problem of connecting the internal fluctuations of a system to its 
response to external forcings and extend the classical theory in order to be able to encompass also nonlinear pro- 
cesses. With this goal, we try to join on the results by Kubo on statistical mechanical systems close to equilibrium, 
i.e. whose unperturbed state can be described by a canonical ensemble, the theory of dispersion relations, and 
the response theory recently developed by Ruelle for non-equilibrium systems equipped with an invariant SRB 
measure. Our derivations highlight the strong link between causality and the possibility of connecting unambigu- 
ously fluctuation and response, both at linear and nonlinear level. We first show in a rather general setting how the 
formalism of the Ruelle response theory can be used to derive in a novel way Kramers-Kronig relations connecting 
the real and imaginary part of the linear and nonlinear response to external perturbations. We then provide a formal 
extension at each order of nonlinearity of the fluctuation-dissipation theorem (FDT) for general systems possess- 
ing a smooth invariant measure. Finally, we focus on the physically relevant case of systems close to equilibrium, 
for which we present explicit fluctuation-dissipation relations linking the susceptibility describing the n th order 
response of the system with the expectation value of suitably defined correlations of n + 1 observables taken in the 
equilibrium ensemble. While the FDT has an especially compact structure in the linear case, in the nonlinear case 
joining the statistical properties of the fluctuations of the system to its response to external perturbations requires 
lineal - changes of variables, simple algebraic sums and multiplications, and a multiple convolution integral. These 
operations, albeit cumbersome, can be easily implemented numerically. 



1 Introduction 

The equilibrium statistical mechanical description of a many-particle system is rooted on the Hamiltonian equations 
of motion, which feature no dissipation, i.e. a vanishing average phase space contraction rate. This basic tenet of 
the Hamiltonian dynamics is mirrored, in the Gibbs' ensemble approach, by the Liouville Theorem, which claims 
the conservation of the probability measure in the phase space. 

On the other hand, nonequilibrium statistical mechanics deals with the investigation of dissipative dynamical 
systems attaining a steady state. Thus, the mathematical description is richer: beside the many-particle system, 
one needs to suitably take into account the effects of the external field performing work on the system, and of a 
thermal reservoir, which absorbs the heat generated within the system by the action of the external field. One of the 
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few general results, in this setting, is represented by the Fluctuation Relations (FR), which concern some peculiar 
symmetry properties of the underlying microscopic dynamics under time-reversal, with respect to the phase space 
contraction rate (which can be replaced, in stochastic dynamics, by the Entropy production of the system). The FR 
hold for systems arbitrarily far from equilibrium and may be essentially regarded as a large deviation result Il26l l5l. 

When considering weakly dissipative systems, the Fluctuation-Dissipation Theorem (FDT) (8l CD establishes 
a link between the response properties of a system to an external (possibly time-dependent) perturbation and a 
correlation function computed in the equilibrium ensemble. The idea underlying the FDT is simple, yet powerful: 
for lineal - deviations from equilibrium, one may compute, for instance, the viscosity of a fluid or the resistance 
of an electrical wire (which are related to the dissipation produced in the system when an external driving, i.e. a 
shear stress or a voltage, is switched on) without actually applying the external field, but just through observing 
suitably defined correlation properties of the fluctuations of the unperturbed system. In other terms, tit is possible to 
establish a correspondence between the internal fluctuations of a system and its response to weak external forcings. 
The practical as well as conceptual implications of this result are immense. 

The problem of determining the response of the system immersed in a thermal bath to external perturbations 
was generalised by Kubo himself in order to study higher order effects on observable, which become practically 
relevant in the presence of stronger forcings. This line of investigation has since start found relevant applications 
in fields such as optics where by the '60s the revolutionary laser technology allowed for studying intense radiation- 
matter coupling and for generating a complex and fascinating phenomenology with wide-ranging theoretical as 
well as industrial relevance. See the book by Bloembergen [38] for a interesting mix of early appraisal of Kubo's 
work and fresh outlook on the very first studies in nonlinear optics, and the book by Butcher and Cotter 1391 for 
a more recent point of view. For rather obvious reasons, optics has historically been the scientific context where 
the response theory has been developed extensively focusing on the frequency domain, rather than on the time 
domain, description. This has led to emphasizing the link between the fact that the system obeys causality when 
responding to an external perturbation and the existence of general integral dispersion relations, commonly known 
as Kramers-Kronig relations, linking the real and imaginary part of the frequency-dependent response to radiation, 
and the related sum rules. This has been of great practical relevance for computing and reconstructing the optical 
properties of matter, natural as well as artificial. An extensive account of this line of work can be found in Il22ll37l . 

Along a different route, the investigation of the response to external perturbations has been extended to en- 
compass the case where the unperturbed system is not at equilibrium, but is rather forced and dissipative, and lives 
in a non-equilibrium steady state (NESS) ll28l . In this state, the phase space continuously contracts, entropy is 
generated, and SRB measures 11161 l30l [T71 provide the natural mathematical framework for describing its statis- 
tical properties. Ruelle [2] recently derived explicit formulas for describing the smooth dependence of the SRB 
measure to small perturbations of the flow in the case of Axiom A systems ifTTll . Such response theory boils down 
to a Kubo-like perturbative expression connecting the terms describing the linear and nonlinear response of the 
system as expectation values of observables on the unperturbed SRB measure. This approach is especially useful 
for studying the impact of changes in the internal parameters of a system or of small modulations to the external 
forcing, and various studies have highlighted the practical relevance of Ruelle theory for studying what we may 
call the sensitivity of the system to small perturbations. In some cases, the emphasis has been on providing con- 
vincing ways to compute the linear response from the unperturbed motion ll36l . in other studies, the authors have 
highlighted the properties of the frequency dependent linear response of the system |29l[33l . Finally, some efforts 
have been directed at extending the analysis of the frequency dependent response to the nonlinear case lPT5l and 
on testing the robustness of the theory with simple chaotic models (31"]. Recently the response theory has been 
used to study the impact of stochastic perturbations P4l and derive rigorously parametrizations for reducing the 
complexity of multiscale systems 1351 . which provides potentially interesting links to the Mori-Zwanzig projection 
theory ll20l . 

The link between linear response of the system to external perturbations and its internal fluctuations is more 
elusive when the unperturbed state is a NESS. In Refs. [2] is is shown that since the invariant measure is 
singular, the response of the system contains two contributions, such that the fust may be expressed in terms of 
a correlation function evaluated with respect to the unperturbed dynamics along the space tangent to the attractor 
and represents the dissipative version of the equilibrium correlation function occurring in Kubo's theory lfT8ll33l . 
On the other hand, the second term, which has no equilibrium counterpart, depends on the dynamics along the 
stable manifold, and, hence, it may not be determined from the unperturbed dynamics and is also quite difficult 
to compute numerically. This expresses the basic fact, already suggested heuristically by Lorenz |[T4l . that in 
the case of non-equilibrium systems internal and forced fluctuations of the system are not equivalent. When 
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devising algorithms for computing the response of the system, in fact, Majda and collaborators are forced to use 
different methods for computing the correlation-like and the additional term described above IT361 . The relevance 
of the novel term spoiling the canonical structure of the FDT for dissipative chaotic systems is still a matter of 
ongoing research. While Majda and collaborators find this term to be of comparable size as the one coming 
from the usual correlation integral, in Ref. [|4l it is shown, by means of low-dimensional solvable models, that 
the novel term introduced by Ruelle is expected to attain its own relevance only in very peculiar situations, such 
as systems with carefully oriented manifolds in phase space and for initial perturbations chosen along the stable 
directions. These observations stem from the fact that physics is mainly concerned with smooth observables 
and with projections from high-dimensional spaces to lower dimensional ones (6l 13. This may explain why 
some attempts of reconstructing the response to perturbations of a complex system such as the climate via the 
application of the classic FDT have enjoyed a good success |[23ll24ll25l . even if it is clear that the performance 
depends critically on the choice of the observable of interest. Moreover, it is important to underline that recent 
works fl2lTl have emphasized that FDT applies for all systems whose invariant measure is smooth, which is, in 
particular, the case for deterministic systems perturbed with noise ll32l . 

The purpose of this paper is twofold: we emphasize the intimate link existing between the causality of the 
response of the system and the possibility of connecting response and fluctuations. First, we wish to show how the 
Ruelle response theory can be used to derive straightforwardly Kramers-Kronig relations connecting at all orders 
of nonlinearity the real and imaginary part of the susceptibility - frequency dependent response of the system 
to perturbations -, and how the susceptibility can be written in terms of unperturbed properties of the system. 
This is accomplished in Sec. |2 Subsequently, we focus on extending the FDT. In Sec. [3] we derive a nonlinear 
generalization of the classical FDT by considering higher orders in the standard perturbative expansion around the 
unperturbed measure, assumed to be absolutely continuous with respect to Lebesgue. In Sec. |4]we present explicit 
calculations for the canonical equilibrium reference frame, which provide the natural extension to arbitrary order 
of the classical Kubo's FDT. In Sec. [5] we present our conclusions and perspectives for future works. Finally, as 
a side note, in App. |A] we show that, when considering perturbations to a canonical ensemble, at all orders of 
nonlinearity the imaginary part of the susceptibility of the observable conjugated to the external field is intimately 
connected to dissipation. 

2 Response theory and Kramers-Kronig relations 

Ruelle [2] recently derived explicit formulas for describing the smooth dependence of the SRB measure of Ax- 
iom A dynamical systems to small perturbations of the flow. Such response theory boils down to a Kubo-like 
perturbative expression connecting the terms describing the linear and nonlinear response of the system as expec- 
tation values of observables on the unperturbed measure. At order n of nonlinearity, such expectation values can 
be written as n— uple convolution of a causal Green function with the time-delayed perturbative fields, so that at 
every order Kramers-Kronig relations can be written for the Fourier transform of the Green function, the so-called 
susceptibility lfT31 . 

Let's consider a general dynamical system whose evolution equation can be written as x = F(x) and let's 
assume that it possesses an invariant SRB measure p(°>. Ruelle lfl9l [T8l [T5l has shown that if the system is weakly 
perturbed so that its evolution equation can be written as: 

x = F(x) + X(x)T(t) (1) 

where X(x) is a weak time-independent forcing and T(t) is its time modulation, it is possible to write the modifi- 
cation to the expectation value of a general observable A as a perturbative series: 

oo 

p{A) t = J2p (n) (A)u (2) 
i=o 

where p(°> is the unpertubed invariant measure, p^ (A) t with n > 1 represents the contribution due to n th order 
nonlinear processes and can be expressed as a n— uple convolution product: 

/OO /'OO 
dn... dr n G(r u T n )T{t - n)T(t - r„). (3) 
-oo J — oo 
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The integration kernel G^(ti, . . . ,r n ) is the n th order Green function, which can be written as: 

G(n, . . . ,<) = J p(°)(dx)G(r 1 )e(r 2 - n) ... 6(r n - T n _i)An(r n - r„_i)An(r n _i - r n _ 2 )An(r 1 )^(x). 

(4) 

where A(») = X ■ V(») describes the impact of the perturbation field and II(cr) is the unperturbed time evolution 
operator such that U(a)K(x) = K(x(a)). The Green function obeys two fundamental properties 

• its variables are time-ordered: if j > k, Tj > Tfc — >• G^ n )(ri, . . . , r n ) = 0; 

• the function is causal: t\ < — > (n, . . . , r n ) = 0. 

Obviously, in the linear case only the second condition applies. These properties allow rewriting the Green function 
as G(t\, . . . , r n ) = 0(Ti)rr™ =2 0(Tj — Tj-i)R(ri, . . . , r n ), where the Heaviside distribution 0(ti) takes care of 
guaranteeing the causality, the terms of the form @(tj — Tj-i) enforce the time-ordering, while R(t\, . . . , r n ) is the 
Response Function, which contains the information about the microscopic dynamics of the system. By applying the 
Fourier transform to Eq. |3l with Y(cj) = FiY (t)) = dt exp[— iut]Y(t) one obtains the following expression 
ESI: 

/OO f'QO n 

dui\ ... / d Wn x (n) (wi,..., Wn )f( Wl )...fK)^-^ Wj ), (5) 
-oo J —oo 3=1 

where p( n ) (A) {(u) = T(p^(A) t ), f(ojj) = F{T{jj)), and the susceptibility x^H^i, ■ ■ ■ ,w n ) is the n- dimensional 
Fourier transform of G^ n \ri, . . . ,r n ) defined as : 



/oo 
drx . . . dr n exp[-iwin] . . . exp[-iw n r n ]G (n) (n, r 2 , . . . 
-oo 



(6) 



while the term containing the Dirac 5 ensures that the frequency of the output is identical to the sum of the 
input frequencies. Note that we use the same definition for the Fourier transform as in (H, while the sign of the 
frequency variable in the integration is opposite to what used in 071 [T31 . which is more common in the optical 
literature. Defining: 

4 n) (n, . . . ,<) = G^( Tl , ...,r n ) + GW(-n, . . . , -r n ) (7) 
G { 2\t u . . . ,r n ) = G^(n, . . . ,r n ) - GW(-n, . . . , -r n ) (8) 



which are different from zero for r n > r n _i > ... > t\ > and r n < r n _i < ... < T\ < and have opposite 

(n) (n) 

parity with respect to exchange of the sign of all of the variables, being G s even and G A odd with respect to 
this symmetry. We have that: 

/oo 
dn . . . dr n exp[-iwiTi] . . . exp[-i6j n r n ]G^ n) (ri, . . . , r n ) (9) 
-oo 

f°° 

2iQ ; {x (n) (wi, • • • ,0J n )} = / dT 1 ...dr n eyip[-iuJiTi]...exp[-iuj n T n \Gy(T 1 



,t, 



(10) 



Thanks to causality Vj we have that 

G( n )(r 1 ,...,T n )=e(T J )G(">(T 1 ,...,r n ) 



9(^)^(7-1,..., r n ). (11) 
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By applying the Fourier Transform to these identities and using the convolution theorem we obtain: 

X (n) (wi, . . . , u n ) = ± (0( Wi )) * (23fl{x (n) (wi, • • • . . . ,co n )} 
1 

1 

~ 2^ 



l - (e(^)) * (2^{x (n) (^i, ■ • • . • • a)}) 
r (e(^)) * (x (n) (wi,...,w 3 ,...,w„)) 

^ (4) +6{U]) ) * (^{x (n) (^i,---,^,---^ n )} 



= (-^(^) + * Qx {n) (^i,---,^,...^n)) (12) 

where V indicates that the integral must be computed considering the principal part and * indicates the operation 
of convolution product. As the same causality argument given in Eq. [TTJcan be repeated for any time variable tj., 
we have that: 

G (n) (n, • ■ • ,r n ) = U^ =1 Q(T n )G^(n, . . . ,r„) (13) 

n^G^jG^Cn,...,^), (14) 

where j\ , . . . , jj. runs over some or all of the indices 1, . . . , n. When taking the Fourier Transform of the previous 
identities, we obtain: 

= n? =1 (-^p + ^ (wii) ^ * (^-x (n Vi, ,o) (15) 

where * must be intended as multiple convolution product of the variables . . . ,0Jj k . Equations [TBI provide an 
alternative expression of generalized Kramers-Kronig relations for nonlinear susceptibilities, first presented in the 
special case of optical processes in 071 . 

We emphasize that Eqs. [§}■ [jjj provide a fundamental connection between the time-dependent response of 
the system to perturbations and the real and imaginary part of the susceptibility. The Kramers-Kronig relations 
establish the correspondence between the fundamental property of causality in the response with the fact that the 
knowledge of only either the real or the imaginary part of the susceptibility is sufficient to reconstruct the full 
frequency dependent response of the system, both in the linear and in the nonlinear regime. 

2.1 Response Function 

We now take a slightly different way for analyzing the frequency dependent response of the system by focusing 
on the the Fourier transform of the response function R(t\, . . . , T n ). We rewrite the definition of the n th order 
susceptibility as follows: 

/oo 
dn... dr n expf-iwm] . . . exp[-iw n r n ]Q(ri) . . . 9(r n - r n _i)i? (n) (ri, . . . ,r n ) 
-oo 

/oo 2 _ n _ 

dux... da n exp[-iuiai] exp[-iw 2 ^ £>j] • • • exp[-w n ^ Oj]Q{o\) . . . Q(a n )S^ n ' (oi, ... ,a n ) 

3=1 3 

/OO n 71 
dux... da n exp[-iai ^Wj] exp[-ia 2 y~)uj] ■ ■ ■ ex.p[-iu n a n ]Q(ai) . . . Q(a n )S^(ai, . . . ,a n ), 
'OO ■ 1 o 



3=1 3=1 

(16) 



where we have performed the change of variables <j\ = T\ and aj = tj — Tj_i Vj > 2 and we have denned 
SV>(* u ...,v n ) = R^(a u ...,Y^<r j ). Defining A( n )(<7i,...,(7 n ) = 0(<7i) . . . Q(a n )SW(ai, . . . , a n ), we 
obtain that: 

n n 
j=l 3=2 

or, in other terms: 

X W(^ - IAJ - 10, . . . , !/„) = AW^,^, ...,!/„). (18) 
Along the lines of the derivation proposed in the previous section, we obtain: 

A (n) ( Wl , ■ ■ ■ ,U>n) = □?=! ("^ ( j-) + ^ K)) * ^--l^A^Ui, . . . , W n )}) 

= n?=i (-^ ( j-) + ^ («*)) * (^{A^K . . .,«„)}) 

= n - B p (4) + 5 K) ) * (^ a(b,(wi ' • • • ' ' ^) 

= nr =1 (~ l -V (j-^j + 5 K)) * (^ W ( W l, • • • , , CUn)) , (19) 

which highlights the fundamental connection, made possible by causality, between the spectral properties of the 
response function and the corresponding susceptibility, at all orders of nonlinearity. Equation [19] suggests that the 
poles of AW(ui, . . . ,LJ n ) are, at most, those of S^ n \u)i, . . . ,oj n ); moreover, the multiple convolution product 
ensures that all the singularities in the right hand side terms which are not compatible with causality are removed. 

Furthermore, if S^> (a\ , a n ) is even with respect to the change of sign of all variables, so that (uii u> n ) 
is even and real, we have that the real part of A^^oji, . . . , u n ) will be given by the sum of all contributions in 

Eq. [19] including an even number of convolutions between S^(oj\, . . . , ,uj n ) and terms of the form V (j^rp 
whereas the imaginary part, conversely, will result from the remaining terms. Note that each time the convolution 
product is applied, the parity of the function is exchanged. Same result will hold if S^(ai, . . . , a n ) is odd: in 
this case, i times the imaginary part of the A^ n \uji, . . . ,u n ) will be given by the sum of all contributions where 
S( n >(u}i, ... ,,oj n ) (which is odd and purely imaginary) is convolved an even number of times with the factors 

^ ( w - ) ' an( * tne rea l °f ( Wl ' • • • ' Wn ) wn l come from the remaining terms. 

In the linear n = 1 case, (w) = A^ (u) and the results presented in Eqs. [15] and [19] lead us to the classical 
results presented by Kubo in the case of perturbations to Hamiltonian systems immersed in a thermal bath JH. In 
fact, we obtain that if S , W(n) is even, it is equal to G^\n), so that its Fourier Transform S^(cji) is equal to 

2K{x (1) (w)},- Instead, if S^(n) is odd, it is equal to (n), and we have that S^(ui) = 2i3={x (1) (w)}. 
Obviously, the parity properties of (n) depend critically on the unperturbed invariant measure and on the way 
the flow is perturbed, and so on the choice of X(x). When nonlinear processes are considered, the link between 
S( n > (ui , . . . , uj n ) and x (^i > • • • > w„) is indeed less trivial, even if Eqs. [T8lfT9"1 provide an algorithmically feasible 
way to unperturbed properties of the system to its response to external perturbations. 



3 Extending the FDT beyond the linear response: general treat- 
ment 

We now wish to explore how to link the (real or imaginary part) of the susceptibility function at various orders of 
nonlinearity to the Fourier transform of correlations of the system in the unperturbed state. In the linear case, this 
is the fundamental content of the fluctuation-dissipation theorem. In order to pursue this line, following Ruelle Q, 
we must assume that the unperturbed invariant measure p(°) (dx) is absolutely continuous with respect to Lebesgue, 
so that it can be expressed as (dx) = (x)dx. In this case, we can rewrite the linear Green function given in 
Eq. [4]for the case n = 1 as a simple lagged correlation between a function C(x) and the observable at a later time 



A(x{t\)) evolved according to the unperturbed dynamics: 



G(n) = J dxp^(x)Q(r 1 )X(x) ■ Vn(n)i(i) = - J dxe( n )pM(x) V ' {p ^^ x)) U(n)A(x) (20) 
dxe(T l )p ( -°\x)C(x)A(x(r 1 )) = [ dxQin)^ {x)C{x{-n)A{x). (21) 



because / dx{V ■ (p ( -°\x)Q(ti)X(x) ■ U(ti)A(x))) = 0, and where we have used the time invariance of the 
measure ~p(°\x)dx in the last step of the derivation. Note that Eq. |2T] provides a very general form of linear 
fluctuation-dissipation theorem for dynamical systems endowed with a smooth invariant measure. 

In order to generalize this procedure for the n th order Green function, we define the adjoint operators for the 
operators A and H(a): 

(a{x),A/3(x)) = {A + a(x),P(x)) 
(a(x), U(a)P(x)) = (U(a)+a(x),f3(x)} (22) 

where the scalar product (•, •) is the ordinary integral evaluated on the support of (x): 



(a(x),P(x))= dxa(x)p(x) (23) 

■V°)(aO>0 

Assuming that such support is compact or that the functions we consider vanish sufficiently fast at infinity, we 
obtain that: 

A(/3(x)) = X(x) ■ V/3(x) -> A + (a(x)) = -V • (X(x)a{x)), (24) 

while the adjoint of the unperturbed evolution operator n(ri) is given by: 

n(cr)(/?(x)) = /3(x(<r)) -)■ U(a) + {a(x)) = a(x(-a)). (25) 

With these definitions, we derive formally from Eq. |4]the following expression: 

fl (n) (n> • • • ,T n ) = (p {0 \x),AU(T n - r„_i) . . . AU( Tl )A(x)) 

= (n(r 1 ) + A + ...n(r n -T n _ 1 ) + A+p(°)(x), J 4( 2; )), (26) 

which gives the n th order response function (and consequently, the Green function) as a n- times correlation. 
Interestingly, when considering Eq. [26l one notes that the dual function II(ti) + A + . . . IT(r n — r n _i) + A + ^°^(x) 
generates the Green functions corresponding to the perturbation flow X(x) (given the unperturbed variant measure 
(x)) for any considered observable,^4(x). Furthermore, following from the definition given in Eq. [T6l we obtain 
the following expression for S^ n \ai, . . . ,a n ): 

S^(ax, ...,a n ) = (pM(x), AU(a n ) . . . AU(ai)A(x)) 

= (Il(a 1 ) + A+ ...IL(a n )+A+p(°\x),A(x)), (27) 

Combining Eqs. [H-Q31 or, alternatively, Eqs. [T6lfT9l with the previous Eq. [26]and considering the content of Eq. 
[49] we obtain a generalized version of the FDT at all order of nonlinear and for rather general statistical dynamical 
systems, specifically for those possessing a smooth invariant measure. In the next section we will show how to 
derive an actual explicit expression for the FDT in the special, albeit most relevant, case of the canonical ensemble. 



4 Extending the FDT beyond the linear response: Canonical En- 
semble 

Following Kubo (U, we now address explicitly the case of an interacting many particle system whose unperturbed 
state is described by the canonical ensemble generated by the Hamiltonian Hq{x), which takes into account only 
the internal degrees of freedom, and analyze the impact of a adding a weak perturbation Hamiltonian H'(x, t) = 



1 



B(x)T(t), where B(x) is an observable conjugated to the external field T(t) fT2"l 191141. In the Kubo framework, 
the perturbed equations of motions can be written as x = F(x) + X(x)T(t), where F(x) = S ■ VHq(x) and 
X(x) = S ■ VB(x), where S is the simplectic matrix. Following the approach highlighted in Ref. |[37l , one may 
adopt a perturbative technique to solve the Liouville Equation for the probability density p t . To this aim, one may 
formally write 

oo 

p« = E^ ■ 

k=0 

This leads to the equation, valid for arbitrary order n: 



dt 



[H ,p^] + [B(x),p^]T(t) (28) 



where, in the classical case [•, •] indicate the Poisson brackets p^ = e denoting the (time-independent) 



equilibrium canonical density, with Z the canonical partition function, and where Eq. (1281) is supplemented with 
the initial condition p t= o = pr^y In the case of a quantum system, we can interpret [•, •] as l/(ih){», •}, where 
{•,•} is the canonical commutator, and p( ) = J2 a l/^exp[— (3E a ]\a)(a\ where the ja)'s constitute a complete 
set of eigenvectors of Hq. As a result, one obtains that the expectation value of a given observable A can be written 
as [271: 

oo 

(A) t = (A) + J2(A) ( t n) , (29) 

n=l 

where we revert to the notation (»)o = Tr{po»}, which is more common in the statistical physical literature, for 
indicating the expectation value po(»), in both the classical and quantum cases. The following expression holds 
for the terms n > 1: 

(Ap = (-l) n [°° dn... dr^TOe^-n) . . . 9(r n -r n _ 1 ) <[J3(-t„), • • • [B{—t\),A] . . .])„ T(t- n ) . . . T(t-r n ) 

(30) 

so that, following Eq. [3]|4l we can express the n th order Green function as: 

G^ B (n, ...,r n ) = (-l) n 9(ri)e(r 2 - n ) . . . 6(r n - r n _i) ([B(-r n ), . . . [B(- n ),A] . . .]) (31) 
where the lower index of the Green function refers to the fact that we are considering the perturbation to the 

(n) 

observable A due to the coupling with the B field, whereas the response function R A B is: 

R%{r u . . . , r n ) = (-1)" ([B(-r n ), . . . [B(- n ),A] . . .]) (32) 
while the variable-wise rearranged function S 

S^ B (n, ...,r n ) = (-1)" ^ [B(- Ty) , . . . [B(—ti),A] ...}) (33) 

' o 

" (n) 

In the following, we will find a compact expression for its Fourier transform S A B , which, combined with what 
discussed in the previous sections, provides the generalization of the FDT in the case of perturbed Hamiltonian 
systems in contact with a thermostat at an inverse temperature /3. 

4.1 Equilibrium correlation functions: Linear case 

Let us now examine, order by order, how a general expression for the FDT emerges from the previous results. We 
first consider the linear response: 

S%(-h) =R%(- h )= -([B(- tl ),A(0)}) 

-i((BH)A(0)) -{A(0)BH)) ) 

~{C A ,B{-ti)-C A , B {-ti-T)) (34) 
in 





Figure 1: Plot of the response function S^' B (t). Shown are the two cases when the latter is odd (red curve) and 



even (blue curve) under time reversal, which entails, respectively, S^ B 

2»{x2i)(w)}. 



(i) 



with r = ih/3, and where CA,B(ti) = Tr{f?(ii)j4(0)/j(°)} denotes the two-time equilibrium correlation function. 
Moreover, in order to obtain the last equality in d34l ). we employed the invariance of the trace under cyclic permu- 
tations and the fact that the operator e^ H ° effects a time translation by the imaginary time (— r). Then, following 
Kubo [8], when going to the Fourier space, it proves convenient to evaluate the complex conjugates of the various 
Fourier transforms. Thus, for instance, one considers: 

\T(C A . B (-t x ))X = [CU,b(-wi)]* = 6a,b{wx) (35) 

where Cab{u\) is the spectral density ifTTI and where the last equality follows from the the fact that Cab(— h) is 
real. By Fourier transforming and by taking the complex conjugate on both sides of Eq. d34l ). we obtain: 

= -t4(1 " e-e**)d AB (ui) (36) 

If, S^ B (ti) is odd under time reversal, as it is commonly assumed, we have, as discussed in Sect. 2, that 

^(wi) = 2i3{x2, ) B )(wi)}, so that: 

&<X%)(toi)} = ^(1 " e-^)C AB ^i)- (37) 

This represents the standard FDT at the first order in the perturbation. It connects the dissipation, related to 
3{x|}g)(wi)}, to the equilibrium fluctuations, contained in Ca,b(wi). On the other hand, if S^ B (tx) is even 

under time reversal, it follows that S^ B (u}i ) = 2K{x? s )(^i)},cf.Fig.[I] 

Starting from the general expression for the response function given in d33l ), which captures the physics of the 
problem, the strategy we employed at the linear order may be straightforwardly repeated at an arbitrary higher 
order, thus leading to a generalization of the FDT for deterministic systems. For practical matters our expansion 
including second or third order around equilibrium is already new and relevant. 



Q 



4.2 Equilibrium correlation functions: second order 



We now present the first extension of the FDT beyond the linear order, thus showing the explict calculations for 
the second order quantities. We start by writing out explictly the expression for the response function: 

4 2 bHi» -*2) = 4 2 b(-*1» "*1 " *2) = (7^2 " *2), A(0)]]> 

1 ;((B(-h - t 2 )B(-t 1 )A(0)) - (B(-h - t 2 )A(0)B(-t 1 )) + 



B(-ti)A(0)B(-<i - * 2 )> + (A(0)B(-ti)B(-ti - t 2 )> ) 

— 2 (Ca,b(— *i) — *i - h) - Ca,b(— *i - t, -ti - t 2 ) - Ca,b(— h, —t\ —t 2 — t) + C A ,b(— *i - t, —t\ - t 2 -(3"$) 



m 

where we have defined the correlation function Ca,b(£i) t 2 ) = (-S(*2)-S(*i)-<4(0))o- Then, we may consider, again, 
the complex conjugate of the Fourier Transform of the three-time correlation functions occurring in Eq. (I38T ). and 
find: 

"+00 



/+00 r+00 
dt 2 / e-^+^C^Hi, -*i - * 2 )rf*i 
-00 J — oo 

/+oo r+oo 
-00 J —00 

with £1 = — ti and £ 2 = — *i — ^2- Thus, one obtains, from Eq.(|38l): 

1 

(-m) 



^(wi, w 2 ) = ^(1 - e-^)(C A , B (wi - w 3 , w 2 ) - e-^- w »)C A ,B(w 2 ,a; 1 - w 2 )) 



= " e"^)^^!,^) (39) 

where we have defined the generalized spectral density: 

cf B {u u u 2 ) := Ca,b(u7i - ^ 2 ,w 2 ) - e-^-^JCA.BC^.wi - w 2 ) . (40) 

Let us now exemplify how we can use these result to reconstruct the response of the system starting from observing 
its fluctuations. We start by considering the observable B(—t\ — t 2 )B{— ti)^4(0) and derive its expectation value 
CA,B(-h, —t\ - t 2 ) = Tr{B{-ti - t 2 )B(-ti)A{0)p^}. We then compute the complex conjugate of the 
two-dimensional Fourier transform of C^b^i,^) and obtain Ca,b(^i — w 2 ,oj 2 ). Using Eq. |40l we construct 
C A B (ui,Ld 2 ) and, eventally, using Eq. [39l we obtain S A B (ui,u 2 ). Furthermore, we plug S ab (oj\,uj 2 ) into 
Eq. [19] and derive, via a double convolution integral, the quantity A A B (ui,u 2 ). Finally, using the definition 

given in Eq. [T8J we eventually obtain %v b( w i > L0 ^)' which contains the complete information on the second order 
response of the system. Therefore, joining the statistical properties of the fluctuations of the system to its response 
to external perturbations requires linear changes of variables, simple algebraic sums and multiplications, and a 
multiple convolution integral. These operations, albeit cumbersome, can be easily implemented numerically. 

4.3 Equilibrium correlation functions: third order 

We hereby present the explict calculations for the third order quantities. As easily seen, the number of terms 
becomes almost unmanageable, but in the next subsection we propose a general formula. We have: 

^a^b(~ ~~ * 2 > — = B-a!b(~ ~~ *i ~~ *2, — *i — t 2 — ts) = 
= -tL([BH ~t2- t 3 ), [B(-h - t 2 ), [B(-h), A(0))]]) 

= — , ..so (Ca, b{— h, —t\ — t 2 , —t\ —t 2 — t 3 ) — Ca,b{— h — t 2 , —ti —t 2 — ts, —t\ — t) + 

- CA,B(~h, —h -t 2 - t 3 , —t\ —t 2 —r) + Ca,b(— h — t 2 — t 3 , —t\ - t 2 - r, —t\ — t) + 

- C A) B(-ii, ~t\ ~ t 2 , — *i - h - t 3 - r) + C A) b(-*i - h, —t\ - t 2 - ts - t, —ti —t) + 

+ Ca,b(-*i, -h - t 2 - t 3 - r, -ti - t 2 - t) - CA,B(-h -t 2 -t 3 - r, -h -t 2 -T, -t\ - t)) (41) 

in 



where we have defined the four-time correlation function Ca,.b(*i>*2>*3) = Tr{B{t 3 )B(t 2 )B(ti)A{Q)p^>}. Sim- 
ilarly to what obtained for the second order terms, we have that: 



/ / e- i ^ h+UJ ^C A ,B(-h,-t 1 -t2,-t 1 -t2-t 3 )dt 1 dt 2 dt 3 (42) 

-oo J —oo J —oo 

Thus, from Eq. (I4T1) one finally obtains: 



S® B (u> u u> 2 ,u, 3 ) = = -^3(1 - e-^ 3 )(C A ,s(wi - w 2 ,w 2 - w 3 ,w 3 ) - e-^ Wl -^)C A ,B(w2- wa.wa.W! - wa) + 
e -^-3) Cjl>B ( Wl _ ^0,3,^ _ W3 ) + e-^ Wl -^)e-^ Wa - W8 )C A , B (w3,W2 -w 3) o;i -wa)) 

1 (1-^3)^(3) (wi)W2jW3) (43) 



(-ihj* y ' A,B 

where we have defined the function 

Cz' B (u)i,u)2,U3) Ca,b{ui -uj 2 ,u 2 -uj 3 ,uj 3 ) - e -W u i- u s)Ci B (aj 2 - W 3 ,W 3 ,Wi - uj 2 ) + 



e 



4.4 Equilibrium correlation functions: general formula 

The results obtained at the lower orders of the expansion pave the way for a straightforward generalization of the 
explicit formulae detailed above. To this aim, let us define the n-time correlation function CA,B(ti,t 2 , . . . , t n ) = 
Tr{B(t n ) . . . B(t 2 )B(ti)A(0)p^}. Then, by induction, it is possible to prove that, for an arbitrary order n, the 
function C^ B {u)\, ...,w n ) attains the structure: 

n— 1 n— 1 jm— 1 32 — 1 

m=0 jm=mjm-l=l ji=l 

(44) 

with 

- = Ur%i, for £ [1, n); 
Wfc 1 cjfc, for k = n. 

For any m € [0, re — 1], the sums on the r.h.s. of Eq.(l44l) yield m ^Zm-i)\ terms (corresponding to all possible 
combinations of time-ordered equilibrium correlation functions), and it is intended, in our notation, that the term 
corresponding torn = yields Ca,b(&i, ■—,&n)- The function C^ b (oji, w n ) generalizes, to an arbitrary order, 

the spectral density Ca,b(w\) appearing in Eq.(l37T). In fact, the function S^ B (ui, ...,ui n ) is easily derived from 
Eq. (|44l and reads: 

S^ B («i , • • • , a*0 = ^ (1 - e"^ )<?W («!,..., «„) (45) 

As described in the subsection dedicated to the second order, it is possible to define an experimental procedure for 
deducing the response of the system from the spectral properties of the observable CA,B(t±,t 2 , . . . , t n ) thorugh a 
cumbersome yet straightforward set of operations. 



5 Conclusions 

The FDT represents a milestone in the endeavour towards a comprehensive theory aimed at connecting the inter- 
nal fluctuations of a system to its response to external forcings. The seminal formulation proposed by Kubo (H 
addressed linear deviations from equilibrium, and has found a vast range of applications in many fields of natural 
sciences. Recent investigations (2J [3] [36] have tried to extend the FDT outside of the original domain by focus- 
ing on linear deviations from nonequilibrium steady states in the framework of of Axiom A dynamical systems. It 

1 1 



has been underlined that when the unperturbed invariant measure is singular with respect to Lebesgue, as usually 
the case in deterministic dissipative systems, a novel term arises, as shown by Ruelle |2j, which seem to prevent 
immediate extension of the FDT around a non equilibrium steady state. This is related to the fact that for these 
systems forced and free fluctuations are not exactly equivalent lfT8ll33ll . yet recent works suggest the possibility 
of practically extending the range of validity of FDT also in this case @ Nonetheless, a response theory can be 
formulated also in this case lfT9l l2l[T8l. and can be successfully framed in terms of frequency dependent suscepti- 
bilities at all order of nonlinearity lfl5l . In this work, instead, we partly pursue a different approach, with the goal of 
highlighting the strong link between causality and the possibility of connecting unambiguously fluctuation and re- 
sponse, both at linear and nonlinear level. We first show in a rather general setting how the formalism of the Ruelle 
response theory can be used to derive in a novel way Kramers-Kronig relations connecting the real and imagi- 
nary part of the the Fourier transforms of the n-th order Green function, i.e. the susceptibility x^H^i, ...,ui n ). 
Moreover, we extend the Kramers-Kronig theory by showing that the application of multiple convolution integrals 
allow to derive the susceptibility from the Fourier transform of corresponding response function S^ n \u}i, ...,u) n ), 
cf. SecJU In this derivation, we shed light on the role of the causality principle (embodied by the sequence of 9 
functions forming the definition of the Green function, cf. Eq. (fT3T >) and of the time symmetries of the response 
function. Thus, Eqs. (fT71) . and ( fl9T ) represent a first, very general, result. Moreover, in the second part of the 
work, we focus on systems whose invariant measure is absolutely continuous with respect to Lebesgue and write a 
formal extension of the FDT to all orders of nonlinearity. We discuss in detail the case of a (classical or quantum) 
Hamiltonian system perturbed by an external field, described by the operator B(t), from its equilibrium state, given 
in terms of the statistical density operator pW . The idea underlying our approach was preliminarily proposed in 
lTT3l in the context of classical Hamiltonian systems. Then, by resorting on a compact general expression available 
for the nonlinear Green function, we succeed to establish a link between the statistical properties of the equilib- 
rium fluctuations, incorporated in the generalized spectral density C^' b (lji, ...,u n ) and the response, related to 
the function X ( w i> ■ ■ ■ j ^n) at an arbitrary order of nonlinearity. In particular, we provide an exact expression 
for the generalized spectral density, Eq. (|44l) . which, supplemented with Eqs. (fT71) and (fl9l ). allows to establish a 
suitable extension of the FDT for nonlinear processes. While the FDT has an especially compact structure in the 
linear case, in the nonlinear case the derivation of the susceptibility of the system from the observation of suitably 
defined correlations requires linear changes of variables, simple algebraic sums and multiplications, and a multiple 
convolution integral. These operations are lengthy but overall trivial and of easy implementation. Concluding, in 
the Appendix we show the the imaginary part of the susceptibility at all orders is related to the power dissipated 
in the system if we select as observable the physical quantity conjugated to the external field. Our method also 
resembles the approach discussed in Ref. flOll , where a compact stochastic version of a generalized FDT is accom- 
plished by means of a perturbation theory around a reference equilibrium state, equipped with a detailed balance 
dynamics. Further connections between the two perturbation theories, employed, respectively, in the deterministic 
and in the stochastic settings, would be worth investigating. 



A Dissipation and Susceptibility 

It is possible to link the imaginary part of the susceptibility to the time-integrated value of a quantity, which, in 
the case of Hamiltonian systems in contact with a thermal bath and in the special context of linear processes, is 
related to the energy dissipation (and consequent entropy production) of the system as a result of the interaction 
with the perturbative field. Hereby, we wish to extend this result in a more general context. We follow the 
formalism presented by Reichl lfT2l . We then consider the quantity p( n \t) = —T(t) ■ d/dtp\ (A), where P(t) is 
a generalized power absorption of the system associated to the interaction between the time varying field T(t) and 
the conjugated observable A. We can rewrite P(t) as follows: 

1 f d ^ 

p( n \t) = —- — / dujdvduji . . . dujit{uj) exp[— iut] — exp[— iut]x . . . , ui n )T(ui\) . . . T(u n )5(v—'S^ uij) 

(2-Ky J dt z — ^ 

(46) 

j ~ n n 

p( n \t) = - — / duoduji . . . duj\T{uj) exp[— iujt](i Uj) exp[— i ujjt]x^ (<jJi, ■ ■ ■ , u n )T(ui) . . . T{uj n ). 

' ^ J 3=1 3=1 

(47) 
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We integrate P^ n \t) over all times: 

/oo r n n 

dtP^(t) = / ^...^^(-^^O^^^Ox^H^i,---,^)^)---^), (48) 

and derive that the total power absorption at order n comes from coupling the response at frequency i Y^=i w i 
with the incoming field at the opposite frequency. Therefore, the larger the band of the time modulation T(t), 
the easier it will be to find matching conditions on the frequency. Since T(t) is a real function, we have that 
T(uj) = (T(oj))*. Therefore, since P(t) is a real function, we derive that: 

/oo p n n 

dt pM(t) = - . . . ,u n )}f (cji) . ..f(u n )f(-Y,«»j), (49) 

J j=l j=l 

which proves that the at all order of nonlinearity, the imaginary part of the susceptibility describes the power dissi- 
pation of the system as defined by p( n )(t) = —T(t) • d/dtp^ (A). In the special case of impulsive perturbations, 
so that T(t) = T 8(t) ->■ f(co) = T , we have: 

/OO p n 

= - / du x . . .dunC^^Wx^Hui, ■ ■ ■ ,O}r re+1 . (50) 
-oo J - =1 

One must note that in the case of a non-continuum spectrum time modulation T(t) - e.g. when T(t) is constituted 
by 2m frequency components (positive and negative) - contributions to the absorption at n th order will come 
only from the generate terms where the sum of n frequencies chosen, possibly with repetition, among the 2m 
frequencies of T(t) match one of the 2m frequencies of T(t) itself. In particular, in the case of a monochromatic 
input T(t), no absorption will take place, e.g., at all even orders of nonlinearity. 
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